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ABSTRACT 

We discuss approximation techniques for estimating spatially varying 
coefficients and unknown boundary parameters in second order hyperbolic 
systems. Methods for state approximation (cubic splines, tau-Legendre) and 
approximation of function space parameters (interpolatory splines) are 
outlined and numerical findings for use of the resulting schemes in model "1-D 
seismic Inversion” problems are summarized. 
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1. Introduction , We discuss here some of our continuing efforts 
on the development of computational techniques for estimation or 
^‘identification" of parameters in second order hyperbolic systems (e.g., 
the acoustic wave equation) . The parameters of interest include bound- 
ary parameters such as coefficients of elasticity and source terms in 
elastic boundary conditions as well as spatially varying moduli of 
elasticity in the partial differential equation itself. Among the 
features of our approach are the following: We do not require an impulse 
or delta function for the source term; indeed the source term need not 
even be parameterized a priori (although it is in the numerical examples 
presented below). Furthermore, the elastic moduli in the system 
equations can be estimated with or without a priori parameterization 
(i.e., assumption of a specific form or shape class). Our ideas are, 
in principle, applicable to vector systems in appropriately defined 
multidimensional domains. 

We combine results from the theory of dissipative operators, linear 
semigroups, approximation theory (splines, spectral methods), and 
optimization techniques in our attempts to develop theoretically sound 
and numerically attractive procedures. While our long term goals 
include development of efficient methods for use in 2-D and 3-D problems 
such as those that arise in "surface seismic" and "bore hole" inversion 
(we have made progress in this direction) , we illustrate some of our 
ideas here with a scalar "model problem" in 1— D. The system in this 
model problem is given by 

(1) pW\t 1^^’ 0<x<l, t>0 

(2) u (t,0) + k u(t,0) = S(t,q) , 

(3) Uj.(t,l) + k2U^(t,l) = 0, 

(4) u(0,x) = 0, u^(0,x) = 0 . 
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Here p is the medium mass density and E is an elastic modulus in (1) 
while the boundary condition (2) is an elastic surface (x=0) condition 
involving a restoring force parameter as well as a parameter (q) 

dependent source term S. This source term is assumed to be the result 
of a perturbing shock to the medium which occurs on the surface at some 
distance from the point of observation. The medium is initially at rest 
coi^ditions (4)) and it is assumed that no waves are reflected 
at a finite lower boundary (x=l); this is given by the absorbing bound- 
ary condition (3) — (this condition can be obtained formally by factoring 

the equation (1) at x=l and taking = / E(l)/p(l)). We make the 

physically motivated assumptions kj^ < 0, k^ > 0 throughout. 

Our model problem consists of using observations of the system (1) - 

(4) to estimate the parameters q = (p ,E,k^,k 2 ,q) . precisely 

we consider the mathematical problem of minimizing the fit-to-data 
criterion 

(5) J(q) = I |u(t^,x ;q) - y..|^ 

i, j ^ 


over a given admissible parameter set Q. Here we assume we have obser- 
vations y^^ for u(t^,Xj) (the "bore hole" problem) — or for u(t.,0) — 

(the surface seismic problem), where u is the solution to (1) — (4) 
corresponding to q. 

It is convenient for our theoretical discussions to ref oriimiiite tills 
problem in terms of an abstract evolution equation in a Hilbert space. 

To do this we first assume that the system (1) - (4) has been trans- 
formed to a system with homogeneous boundary conditions 

^'^tt “ ?x 

v^(c,0) + k v(t,0) = 0 
(6) "" 1 

v^(t,l) + k 2 V^(t,l) = 0 

v(0,x) = <})(x,q), v^(0,x) = V;(x,q) . 


This can be done in a straightforward fashion (see [ll]). 
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As our state space we choose Z = H (0,1) 
dependent inner product 

1 

<z,w> = / Ez'w'dx - E(0)k. z. (0)w (0) 

’ q 1 1 11 1 


X H^(0,1) with parameter 
1 

+ / pz„w_dx , 

0 


for z = (z^,Z 2 >, w = (w^jW^) in Z. Under reasonable assumptions on 

p, E and k , this yields a Hilbert space Z = Z(q) with topology equiv- 
^ 10 

alent to the usual H x H topology. The transformed system (6) can 
then be written in abstract form as (we don' t distinguish between a 
vector and its transpose) 

z(t) = A(q)z(t) + G(t,q) 

(7) 

z(0) = $(q) 


where z = (v,v^), $ = (<{i,ijj), G = (0,g), and the operator A(q) defined 
on dom(A(q)) = {z € x | zj^(O) + kj^z^(O) = 0, z^d) + k2zj^(l) = 0} 
is given by 


A(q) 


■f 



It can then be shovm that under boundedness assumptions on Q, there 
exists a constant oj independent of q in Q such that A(q) - col is 
dissipative in Z (i.e., <A(q)z,z> < co<z,z>). Furthermore A(q) generates 
a strongly continuous semigroup S(t;q), t ^ 0, that is the family of 
solution operators for (7). 

This framework provides a convenient setting for discussion of 
semidiscrete approximation schemes (and their convergence properties) 
for solving the problem of minimizing (5) over Q. In the case under 
consideration here both the state space for (1) - (4) and the parameter 
set Q are infinite dimensional. We first turn to a description of 
state approximation ideas that we have employed. 
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Abstractly, one approximates (7) — at least in its state variable — 

by choosing a sequence of finite dimensional subspaces Z^, N=l,2,..., 

of the state space Z. Letting P be the projection of Z onto and 
N N N 

A : Z Z be a family of approximating operators for A, one can 
define the sequence of approximating systems in Z^ by 


( 8 ) 


z^(t) = A^(q)z^(t) + P^G(t,q) 
z^(0) = P^<I> 


and the corresponding fit-to-data criterion 



The problem of minimizing over Q is then a finite dimensional state 
problem which can in some cases (e.g., when the functional parameters 

are assumed in parameteric form) be readily solved for approximate 
-N 

parameters q , N = 1,2,*.. . In this situation one then desires to 

argue that the sequence {q } (or some subsequence) converges to a 
parameter q* in Q that provides a minimum for (5). However, in other 
cases where the parameter functions (such as p and E in (1)) are not 
assumed to possess a priori finite dimensional parameterizations , the 
problems involving (8), (9) still entail optimizations over an infinite 
dimensional set and thus a parameter approximation scheme must further 
be introduced before the approximating problems are easily solved. 

We defer discussion of this aspect of our methods until section 4 
below. Instead we first discuss two state approximation schemes 
(cubic spline and tau-Legendre) that we have used with some success. 

2* Cubic spline state approxim ations. These methods are founded 

on ideas given in [l],[3],[9] where cubic B-splines Bj (see [l3] for 
general concepts and discussions) are used to generate basis elements 
N 

for subspaces Z in which the elements satisfy the boundary conditions 
inherent in the problem. For the problems under consideration in this 

note, this results in parameter dependent basis elements = B^(q). 

Thus the dependence of = Z^(q) on q is not only because of the 
parameter dependent inner product but also through the explicit 
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dependence of the basis elements on the parameters. This aspect 

leads to nontrivial technical difficulties (from both a theoretical 
and a numerical viewpoint) in extending the ideas of [ 3 ]. One can 
however overcome these difficulties to obtain a theoretically sound 
and computationally feasible scheme (see [s], [ll] for more detailed 
discussions) • 


To give a brief idea of the theoretical aspects of this method, we 


N N 

first observe that the above considerations lead to subspaces Z = Z (q) 
in the context of the Galerkin approach of [ 3 ], [?]. One then desires 

to define Galerkin approximations = A^(q^) = P^(q^)A(q^)P^(q^) 

N 

associated with a given sequence of parameter estimates {q }. Here 

N N 

P (q) is the orthogonal projection of Z(q) onto Z (q) • To establish 

a state convergence theory for (8), it suffices to argue stability 

and consistency for the schemes (see [3], [?]). Stability follows 

readily from the uniform dissipativeness of the operators 


A(q) : <A (q)z,z> = <P (q)A(q)P (q)z,z> = <A(q)P (q)z,P (q)z> < 

N N 

o)<P (q)z,P (q)z> < 0 )<z,z>. Consistency is somewhat more delicate 

N * 

since one wishes to argue that for q q in Q one has 
N N * 

A (q )z A(q )z. More generally in the process of discussing con- 
sis tency and state convergence one must deal with statements Involving 

N N N 

the convergence of elements z in Z (q ) — which satisfy boundary con- 
, N * * 

ditions involving q — to elements z in dom(A(q )). Hence one must 
construct some device to express convergence of the boundary condition 
parameters as well as that of the parameters defined directly in the 
N N 

operators A (q ) . With care this can be done and one can employ the 
Trotter-Kato approximation theorem from linear semigroup theory [l2] 
along with estimates from spline approximation theory [l3] to establish 
consistency and state convergence (see [S], [ll] for details). 


These spline based schemes have proven quite satisfactory in 
numerical computations with test examples as we shall see in the 
numerical section below. 


3. Tau-Legendre state approximations . The tau method (due to 
Lanczos) is a special case (see 11 of [lO]) of the spectral methods 
discussed in [lO] and involves use of eigenfunctions that are in 
general not related to the natural modes of the system being approxi- 
mated. As outlined above, one again rewrites the system (6) as an 
abstract system in the state space Z = Z(q) defined in section 1. 

The system (7) is written in the form 


( 10 ) 


z(t) = L(q)z(t) -h G(t,q) 
z(0) = $(q) 
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(11) B(q)z(t)=0, 


where L(q) is the same as the operator A(q) in (7) except that we do 
not include the boundary conditions in defining the domain of L(q); 

2 1 

i.e. , dom(L(q)) = H x H . Instead the boundary operations are 
imbedded in an operator B and we impose the side conditions (11). 


The approximating subspaces Z are then defined in terms of an 

orthonormal family that is complete in H^(0,1). Members in this family 
are not required to individually satisfy the boundary conditions; these 
are imposed on the approximations to the solution of (10), (11) by 
requiring that the approximations satisfy exactly the conditions (11) . 
To be more specific, suppose we begin with N+1 elements 


(e.g., Legendre functions) to generate the appropriate subspaces 
N N 

Z (q) of Z(q) — (the Z (q) have dimension 2N = (N+1) + (N-1) in this 
case — the analogous cubic spline approximating subspaces have dimension 
2N+3). Then we assume that the first components of solutions z= (v,v ) 

of (10) — i.e., of (6) — are approximated by v (t,x) = ^ w.(t)((i.(x) 

j=l ^ ^ 

where the first N-1 coefficients w^ are determined by Imposing the 

equation (10) and the remaining coefficients are determined by imposing 
the boundary conditions (11) . 


We observe that the tau method is not a Galerkin procedure. It also 
differs from the cubic spline method described above in that the 

N N+K 

approximating subspaces satisfy Z C Z , K > 0. 

Our efforts have been devoted to the use of the tau-Legendre method 
in the model problem of section 1. In this case, we chose 
(pj (x) = Pj_^(2x - 1) , 0 _< X £ 1, j =1,2,... ,N+1, where is the 

Legendre polynomial of degree j. For a convergence analysis, one can 
employ the same Hilbert space framework described in discussing the 
cubic spline schemes. In this operator formulation one again estab- 
lishes convergence via demonstrating stability and consistency. The 
arguments involve dissipative estimates, the Gronwall inequality, and 
approximation properties of Legendre polynomials. 


One interesting (and potentially highly advantageous) aspect of 
the tau-Legendre methods involves their use in layered media problems 
(e.g., see Example 4 below) where the coefficients p, E in (1) possess 
discontinuities. In this case one can include (in addition to the 
boundary conditions) in the conditions (11) the continuity conditions 


N - N + - N - 

V (t,x. ) =v (t,x. ) — on displacement — and E(x. )v (t,x, ) = 
J J J ^ J 


+ N + 

E(x. )v (t,x. ) — on stress — at interfaces x, in the medium. 

J X J J 


This 
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could lead to methods that are computationally superior to methods in 
which these interface continuity conditions are imposed directly on the 

N 

basis elements for Z . 


4. Parameter approximations . As we have noted above, the problem 
N 

of minimizing J of (9) over Q subject to the state approximatxon 
equations (8) is still a difficult infinite dimensional optimization 
problem in cases where Q has nonparameterized function space components 
(e.g., in situations where we estimate p and/or E in (1) choosing from 
some infinite dimensional function classes). In such situations it is 

M 

useful to introduce a family of finite dimensional parameter sets Q 
that approximate, in some sense, the original parameter set Q — i.e., 

"Q^ ^ Q" as M -»■ «. We may then consider the problem of minimizing 

of (9) over Q^, which is finite dimensional in both the states and 

M 

parameters and may be easily solved for judicious choices of Q , 

, . . -N 

producing estimates . 


We have, with some success, used such double approximation procedures 
to estimate functional coefficients in parabolic transport equations 
[^]. [ 5 ], [6] and in higher order equations of elasticity [2], In 

these efforts we employed approximating sets Q consisting of either 
linear or cubic interpolatory splines. Under reasonable assumptions 

on Q and on the approximation properties of Q (e.g., see [2], [4], 

[6] where theoretical results are also given for the systems mentioned 
above), one can prove a double limit convergence result: Given any 



as defined above, there exists a subsequen- 
Q that provides a minimum for J of (5) 


In the next section we shall present a sample of some of our 
numerical results using such procedures. 


5, Numerical examples . We present here a summary of some of our 
numerical findings for the methods described above when used with the 
model problem involving (1) - (5) . The estimation schemes were tested 
on examples in the following manner. First, known ”true” values 
(denoted by q* in the tables below) of the parameters to be estimated 
were chosen and a numerical method independent of the one being tested 
was used to solve the forward problem for a solution u = u(q ). Values 
of this numerical solution were used for ’’observations” or data for 
the inverse algorithm. The minimization problems were solved itera- 
tively (with "start-up” or initial guess q^ listed in each example) 
employing an IMSL version (ZXSSQ) of the Levenberg-Marquardt algorithm. 
The ordinary differential equations for the approximate states were 
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solved using either DVERK or DGEAR in the IMSL packages. The test 
calculations were carried out on an IBM 370 at Brown University, a 
CDC 6600 at Southern Methodist University, and CDC Cyber 170 at*NASA 
Langley Research Center. In the first two examples we compare the 
cubic spline state approximation scheme with the one based on tau- 
Legendre state approximations. In the third example we give a sample 
of our findings on a problem where the double approximation ideas of 
section A were used while in the last example we detail results for a 
multilayered media problem. 

Example 1 . In this example we considered (1) - (4) with E/p = q 

constant, source S(t,i) = q (1 - e"^*^)e'^^ and k, , k. to be estimated. 

The "true" values q , start-up values q^, along with findings for the 
cubic spline state approximations and the tau-Legendre state schemes 
are given^below in tabular form. Observations included values 

) for t. values .5, land 1.5, x. values 0, .5, and 1. The 
7 '* N -N J 

residual (J (q )) for the converged parameter values along with CPU 

time in seconds is listed for each level (N) of state approximation. 



-N 

‘^l 

-N 

^12 

-N 

‘*3 

k^ 

‘"l 

k^ 

*"2 

jN.-N. 

J (q ) 

CPU 

CUBIC SPLINES 








N = 4 

3.014 

2.001 

-1.084 

-1.883 

.992 

. 6x10“^ 

124 

N = 8 

2.992 

2.017 

- .958 

-2.082 

.999 

.17x10"^ 

108 

TAU-LEGENDRE 








N = 4 

2.940 

2.098 

-1.026 

-2.034 

.994 

.7x10 ^ 

31 

00 

II 

3.017 

1.957 

-1.001 

-1.979 

1.008 

. 6x10”^ 

40 

TRUE ^ 

VALUES (q ) 

3.0 

2.0 

-1.0 

-2.0 

1.0 



START 
UP (q^) 

2.0 

1.5 

- .5 

-1.0 

2.0 




Example 2 . We considered the transformed system (6) with p = 1, 
g = 0, parameterized modulus E(x) =1.5+1 Arctan(q^(x- q^)) and* 

Loitial data <t>(x) = e* , ^ (x) = - 3e^. Observations at values 
f -33, .5, .66, .83, 1 and x^ = 0, .5, 1 were used. The 

following converged values were obtained. 
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-N 

‘’l 

-N 

^2 

k^ 

k^ 

2 

.N,-N. 
d (q ) 

CPU 

CUBIC SPLINES 







N = 4 

2.964 

.487 

-.990 

3.006 

.8x10"“^ 

60 

N = 8 

3.051 

.501 

-.999 

2.996 

.28x10"^ 

54 

N = 16 

3.012 

.5002 

-.9999 

2.999 

.19x10"^ 

146 

TAU LEGENDRE 







N = 4 

3.191 

.516 

-1.028 

3.015 

.12x10"^ 

14 

N = 8 

2.942 

.496 

-1.001 

2.997 

.8x10 ^ 

57 

N = 16 

2.991 

.4992 

- .9998 

2.999 

.15x10"^ 

520 

TRUE ^ 

VALUES (q ) 

START 
UP (q^) 

3.0 

.5 

-1.0 

3.0 



5.0 

1.0 

-2.0 

2.0 




Example 3. We considered the system (6) with p - 1, g - 0, 

(j)(x) = e^, ijj(x) = - 3e^ and parameters to be estimated E (x) = 1.5 + 

* * . • 
tanh(6(x- .5)), kj^ = -1.0, k .2 = 3.0. We did not assume an a prion 

parameterization for E, but used linear interpolatory splines for the 

parameter approximation procedure as described in section 4. Cubic 

spline state approximations were used. Observations were taken at 

the same t., x, as in Example 2. The initial guesses were E^(x) = 1.0 

^ ^ 0 0 -N . 

(a constant function), k^^ = - 2.0, k 2 = 2.0. Graphs for E^^ are given 

in the figures below for N = 4 (2N+3= 11 cubic elements in the basis 

for Z^) and N = 16 with parameter approximations corresponding to 
M=3 (M+l = 4 linear elements to approximate E) . For N = 4, the 

results obtained were k^ = - 1.054, k^ = 3.357, J (q ) = .25x10 , 

[E - E^I = .81x10 , CPU = 38 sec., while for N = 16 the results were 

= - 1.100, k 2 ^ = 3.070, = .47x10"^, | E* - E^^j = .29xl0"^, 

CPU = 118 sec. 
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Example 4 . We considered system (1) - (4) for a multilayered medium 

(i.e., a discontinuous elastic modulus). We chose p = 1, E piecewise 

^ ^ 3 1 — t 

constant (E given below) and S = 2(l-e )e and used the tau- 

Legendre method for the state approximations as described in section 3 

to estimate E, k^, k 2 * Observations consisted of surface observations 

(Xj = 0) at 40 evenly spaced times t^ in [0,2] (i.e., t^ = .05, .1, 

.15, ...). For N = 4, we obtained the estimates listed below with 

J^(q^) = .3x10"^, CPU time = 172 sec. 


INITIAL GUESS 



"1.5 

0 £ X 

< 

.2 


II 

0 r-C 

-2.0 

E°(x) 

1.5 

.2 < X 

< 

.8 





1.5 

.8 < X 

< 

1.0 



2.0 

TRUE VALUES 

^2.0 

0 _< X 

< 

.3 


* 

^1 ^ 

- 1.0 

* 

E (x) = 

3.0 

. 3 < X 

< 

.7 





^1.0 

. 7 < X 

< 

1.0 


■k 

\ = 

3.0 

ESTIMATES (N = 4) 








^1.996 

0 

< 

X £ 

.299 

^1 " 

- .998 

E^(x) = ^ 

2.986 

.299 

< 

X £ 

.7006 




^ .996 

.7006 

< 

X £ 

1.0 

= 

2.987 


Concluding remarks . As we have outlined above, a convergence 
theory can be developed for both the state and parameter approximation 
schemes described in this note, ^e have tested the methods on a 
number of examples in addition to the ones presented here. The schemes 
perform well on both "bore hole" and "surface seismic" type model 
problems. In some cases (see Example 1) the tau-Legendre state 
approximations perform slightly better than the cubic spline approxi- 
mations. In others (see Example 2), the cubic spline schemes appear 
to be more desirable. 


While our test examples are clearly very simple, they do serve the 
purpose of illustrating the potential for use of the methods and ideas 
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in estimating variable (continuous and/or discontinuous) coefficients 
and unknown boundary parameters in second order hyperbolic systems. 
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